Virial expansion for a strongly correlated Fermi gas 
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Using a high temperature virial expansion, we present a controllable study of the thermodynamics 
of strongly correlated Fermi gases near the BEC-BCS crossover region. We propose a practical way 
to determine the expansion coefficients for both harmonically trapped and homogeneous cases, and 
calculate the third order coefficient 63 (T) at finite temperatures T. At resonance, a T-independent 
coefficient 63°™ » -0.29095295 is determined in free space. These results are compared with a recent 
thermodynamic measurement of 6 Li atoms, at temperatures below the degeneracy temperature, and 
with Monte Carlo simulations. 
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Strongly correlated Fermi gases are of wide interest and 
underlie many unanswered problems in quantum many- 
body systems, ranging from neutron stars, hadrons and 
quark matter through to high T c superconductors [l| . Re- 
cent investigations of Feshbach resonances in ultracold 
atomic Fermi gases have opened new, quantitative op- 
portunities to address these challenges [l]. A great deal 
of theoretical work has been carried out for this simple, 
well-controlled case of a strongly interacting yet low den- 
sity Fermi gas, which is known as the unitarity limit. 
However, a profound understanding is plagued by the 
large interaction strength, for which the use of pertur- 
bation theory requires infinite order expansions. Numer- 
ically exact quantum Monte Carlo simulations are also 
less helpful than one might expect @, H, 0]. Due to the 
Fermi sign problem^, computer simulations are often 
restricted to small samples, and are therefore difficult to 
extrapolate to the thermodynamic limit. 

In this Letter, we approach this problem by using a 
controllable virial expansion study of trapped strongly in- 
teracting Fermi gases at high temperatures. We focus on 
the low-density physics which is described by an effective 
S-wave contact potential. Our expansion has a small pa- 
rameter. The fugacity 

z = cxp(/i/fc B T) <C 1 

is small because the chemical potential fi diverges loga- 
rithmically to —00 at large temperatures T. The virial 
expansion up to the second virial coefficient was applied 
by Ho and Mueller to explore the universal thermody- 
namics of a homogeneous Fermi gas at unitarity [H| • Here 
we extend this to the third order coefficient. Most impor- 
tantly, we present a practical theoretical strategy which 
can even be extended beyond third order. Surprisingly, 
we find that the simplest theoretical route to calculating 
these higher order coefficients is via the use of exact solu- 
tions for the energy eigenstates of harmonically trapped 
clusters. This gives a unified approach to calculating 
virial coefficients in both trapped and untrapped cases. 




Figure 1: (color online). Predicted interaction energy as a 
function of entropy at unitarity, compared with experimen- 
tal data from Duke University f 1 01 ] - The solid (dashed) line 
shows the contribution up to the third (second) order virial 
coefficient. The arrow indicates the degeneracy point. 



In reality, the strongly interacting spin- 1/2 6 Li and 
40 K fermionic gases were realized by tuning a magnetic 
field across a resonance []]]. Extensive experiments have 
studied the crossover from the BCS limit (Cooper pair- 
ing of atoms) to the BEC limit (Bose-Einstein conden- 
sation of diatomic dimers). The most interesting region 
lies at the middle of crossover, where the two-body S- 
wave scattering length a becomes much larger in magni- 
tude than the inverse Fermi vector 0. Fascinating 
phenomena may occur in this 'unitarity limit' [3,13], such 
as the observed scale-invariant, universal thermodynamic 
behaviour 0, [13, 11 1. Accurate, high-order virial coeffi- 
cients provide an extremely useful tool in analysing these 
experimental results above the superfluid transition. At 
the same time, the methods given here may have general 
applicability to other strongly interacting systems. 

In what follows, we first introduce a practical way to 
calculate the n-th virial coefficients b n (T) and, by solving 
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exactly the two- and three-particle problems, determine 
the second and third virial coefficients in an isotropic har- 
monic trap at the BCS-BEC crossover. We then focus on 
the unitarity limit and calculate the energy and entropy 
of a trapped gas using the virial expansion method. 

Our main result is summarized in Fig. 1, which shows 
the comparison of the virial expansion prediction to a 
recent measurement of the entropy-dependence of the in- 
teraction energy E int (S) at unitarity. The experiment 
was carried out for atomic gases of 6 Li atoms at a broad 
Feshbach resonance (lp| |. We find an excellent agreement 
at temperatures below the Fermi degeneracy temperature 
Tp. This remarkable result is opposite to the consensus 
that the virial expansion is valid at the classical Boltz- 
mann regime with T ^> 1>. We suggest that it can be 
understood by the significant suppression of higher-order 
virial coefficients in a harmonic trap. 

Even in the absence of a harmonic trap, our method 
can still be used to calculate the virial coefficients. We 
determine a universal coefficient o|°™ ss —0.29095295 for 
a homogeneous Fermi gas at unitarity, in contrast to a 
recent calculation that obtained a result with the oppo- 
site sign [13]. Our resulting equation of state is in good 
agreement with existing Monte Carlo results, and may 
provide a useful benchmark for testing future quantum 
Monte Carlo simulations of strongly interacting Fermi 
systems at high temperatures. 

Virial expansion. — Let us consider the ther- 
modynamic potential = —kBTlnZ, where Z = 
Tr exp[— (7i — /iA/")/fcsT] is the grand partition func- 
tion. At high temperatures, we can rewrite Z in 
terms of the partition functions of clusters, i.e., Q n = 
Tr„[exp(— TLn/kBT)] with n denoting the number of par- 
ticles in the cluster and Tr n denoting the trace over 
n-particle states of the proper symmetry; thus we find 



Z = 1 + zQi + z 2 Q 2 + 
tial can then be written as, 



The thermodynamic poten- 



-k B TQ\ 



\z + b 2 z 2 + 



+ b n z n + 



where the virial coefficients are given by, 

b 2 = (Q 2 -Q?/2)/Qi, 

h = (Q 3 - Q1Q2 + <2i/3) /Qi, etc. 



(1) 



(2) 
(3) 



These equations present a general definition of virial 
expansion and is applicable to both homogeneous and 
trapped systems. The determination of the n-th virial 
coefficient thus requires full solutions up to the n-body 
problem. It is convenient to focus on the interaction ef- 
fects only and consider Ab n = b n — bn and AQ n = 
Qn — Qn , where the superscript "1" denotes the non- 
interacting systems. We shall calculate A6 2 = AQ 2 jQ\ 
and Ab 3 = AQ3/Q1- AQ 2 . 

Second and third virial coefficients in traps. — By 
solving the few-body problem exactly, we now evaluate 



the virial coefficients in a three-dimensional isotropic har- 
monic potential V(r) = miur 2 /2, with a trapping fre- 
quency lj and fermion mass m. The partition function 
Qi is easily obtained from the single-particle spectrum of 
the harmonic potential, E n i = (2n + I + 3/2)Hlj, and the 
single-particle wave function, R n i{r)Y l m {6, ip). We find 
that Qi = 2exp(— 3w/2)/[l — cxp(— u))] 3 with a dimen- 
sionless frequency u> = huj/ksT -C 1. The prefactor of 2 
in Q\ accounts for the two possible spins of each fermion. 

To solve the two- and three-fermion problems, we 
adopt a short-range S-wave pseudopotential for inter- 
atomic interactions, in accord with the experimental situ- 
ation of broad Feshbach resonances. This can be replaced 
by the Bethe-Peierls contact conditions on the wave func- 
tion tp (ri, T2, r n ); when any particles i and j with un- 
like spins close to each other, ry = |r^ — — * 0, r^tp 
satisfies, 



d{rijip) /drij = - (njip) /a. 



(4) 



Otherwise, the wave function ip obeys the noninteracting 
Schrodinger equation, 



h 2 , 

2m r 



1 ! 

—mujr 4 
2 



tp = Etp. 



(5) 



Unlike Bose gases, no additional many-particle interac- 
tion parameters are required. This is due to the absence 
of many-particle S-wave bound states (i.e., Efimov-like 
states) for these low density Fermi gases, which has a 
physical origin in the Pauli exclusion principle. 

The Hamiltonian of two fermions with different spins 
was solved by Busch et al. As the center of mass 

is separable for a harmonic trap, we may single out the 
(single-particle) center-of-mass energy £7 c . m , and rewrite 
the total energy as E = E c . m . +E re i- Following Busch et 
al. the relative energy E re i — (2v + 3/ 2)fko satisfies, 
2T(-v)/T(-v - 1/2) = d/a, where d = ^2h/muj is the 
length scale of the trap, and the (un-normalized) two- 
body relative wave function is given by, tp^f (r = rj — 
rr,v) = exp(-r 2 /2rf 2 )r(-z/)J7(-^,3/2,r 2 /d 2 ). Here, 
the total angular momentum of tp^f is strictly zero since 
only these states do not vanish at r — and thus 
are influenced by the pseudopotential. T and U are 
the Gamma function and confluent hypergeometric func- 
tion, respectively. It is readily shown that 62 — fr 2 = 
AQ2/Q1 is given by: 6 2 -6 2 1} = (1/2) £^[ e -(2-n+3/2)a_ 
e -( 2ly i ) +3/2)o)j^ wnere the- summation over £ c , m , cancels 
Qi in the denominator, and v n ^ = 0,1,... is the n- 
th solution of the relative energy spectra in the non- 
interacting limit. At unitarity, the two-body solutions, 
= n — 1/2, are known exactly [13], leading to 
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2 [1 + exp (-£)] 4 32 



(6) 



The three-fermion problem was studied analytically by 
Werner and Castin [l4j at unitarity, and numerically by 
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Figure 2: (color online). Relative energy levels of a three- 
fermion system at the ground state section (1 = 1). 
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Figure 3: (color online). The second and third virial co- 
efficients as a function of the interaction parameter l/k F a. 
We have used a total number of atoms N = 100, leading to 
Q = (3AT 1/3 w 0.15 at T = T F . 



Kestner and Duan [15] for arbitrary scattering lengths. 
Although the calculations are more involved, the exact 
solution is intuitively understandable. Let us skip the 
trivial center-of-mass motion. Assuming a spin state it! 
and using the Jacobi coordinates r = r 2 — ri and p — 
(2/v3)[r3 — (ri/2 + ra/2)] as shown in the inset of Fig. 
2, the three-body relative wave function ip^ (r, p) can be 
written as 

r 3 f = (1 - V 13 ) J2 a n R nl (p) Yr (p) V 2 t (r; u n ), (7) 

which is simply the summation of products of the 
eigenstate of the paired particles 1 and 2, i/ , 2fe'( r ! u n), 
and of the eigenstate of particle 3 relative to the 
pair, R n i(p)Y [ m (p). The value of u n for each index 
"n" is uniquely determined from energy conservation: 
Erei/ha; = (2n + l + 3/2) + (2u n + 3/2) and should not 
be confused with the solutions for the two-body rela- 
tive energy. The relative wave function tp^f has a well- 
defined total relative angular momentum of the 3 parti- 
cles with quantum numbers I and m. The operator Vn 
ensures the correct exchange symmetry of the wave func- 
tion. This introduces correlations between the 1 — 2 pair 
and the remaining particle 3 and thus a hybridization as 
parametrized by a n . We solve the eigenstate "a n " and 
eigenvalue E re i by imposing the Bethe-Peierls boundary 
condition Eq. |(4]). We find, 



2r(-o 



1/2) 



Cn 



- | a n , 

a, 



(8) 



where the (symmetric) matrix C nm = [(— 1)' / 'V^\ x 
J °° dpp 2 R ni (p) R ml (p/2) ipl z f(v / 3p/2; v m ) arises from 
the exchange operator P13. Without C nm we have a 
three-body problem of un-correlated pair and single par- 
ticle. We label the relative energy in this case as Erei and 
calculate it directly from the two-body relative energy. 

We have solved Eq. ([8]) numerically for 10 4 energy 
levels E re i at different relative angular momenta I and 



have checked that at unitarity our results agree exactly 
with the analytic spectrum in Ref. [1J|, with relative 
numerical errors typically < 10 -6 . Fig. 2 shows how the 
relative energy spectrum evolves from the BCS to the 
BEC side in the subspace of I — 1. 

To calculate the third virial coefficient using 63 — 63 = 
AQ3/Q1 — AQ2, we notice that the spin states of if J, 
and tit contribute equally to Q3. Also, Q\ in the de- 
nominator is canceled by the summation over E c , m ,, and 
the term — AQ2 is cancelled by the difference between 
E re i and the noninteracting energy E^J. Thus, the third 
virial coefficient is determined solely by the exchange 
correlation, so that 63 — = ex P (~ E re i/kBT) — 
J^exp {—Erei/hsT), where the summation is performed 
over all possible three-body states that are affected by 
interactions. At unitarity, we obtain: 



^3,00 



y 3,oo 



-0.06833960 



0.038867w 2 



(9) 



The second and third virial coefficients through the 
crossover is given in Fig. 3 at three typical tempera- 
tures. Here we consider a gas with N = 100 atoms and 
scale the inverse scattering length using the Fermi vector 
at the trap center, k F = (24A^) 1 / 6 /(d/ % /2). The tem- 
perature is given in units of Fermi temperature T F = 
E F /k B = (3N) l / 3 (hu>/k B ). All the curves with distinct 
temperatures cross at a — * ±00. This is the manifesta- 
tion of universal behavior anticipated if there is no any 
intrinsic length scale. However, the characteristic length 
scale d of harmonic traps brings a small (non-universal) 
temperature dependence that decreases as N~ 2 / 3 , shown 
by the terms u> 2 in Eqs. ((6]) and 

High-T thermodynamics in traps. — We are ready to 
investigate the thermodynamics of a strongly interacting 
Fermi gas at high temperatures. At unitarity, the energy 
E and entropy S in the limit of small Q or large N can 



4 



be calculated according to the universal relations @, S], 



E = -3aQ/2, 
S = -(3a/2- 



l)n/T-k B Nhiz, 



(10) 
(11) 



together with Eq. |T]) for fl and the number identity 

N = - (dfl/dfx) = Qi [z + 2b 2 z 2 H ] , and the virial 

coefficients shown in Eqs. J6|) and ([9]). Here, a = 2 for a 
harmonically trapped gas. We also calculate the equation 
of state Eig(S) of an ideal non- interacting Fermi gas, 
using b { n ] = [1/n 4 - £ 2 /(8n 2 ) + ■••]. 

Fig. 1 shows the predicted interaction energy Ei nt = 
E — Eig as a function of entropy, compared to the ex- 
perimental data reported by Luo et al. [IB]. We find a 
rapid convergence of expansion, even below the degener- 
acy temperature Tp, with excellent agreement between 
theory and experiment: the virial expansion is applicable 
to a trapped Fermi gas even at T <Tp. 

This remarkable observation is counter-intuitive, as 
the virial expansion is generally believed to be useful 
at the Boltzmann regime with T 3> Tp. This oc- 
curs because there is a significant reduction of higher- 
order virial coefficients in harmonic traps. Consider the 
thermodynamic potential of a harmonically trapped gas 
in the local density approximation, f2 = J il (r) dr oc 
/ dr[z (r) + b*°™z 2 (r) + ■■■ + b^™z n (r) + •••], where 
z(r) = zexp[— V (r) /ksT] is a local fugacity with the lo- 
cal chemical potential /z(r) = fi— V (r). It is readily seen 
on spatial integration that the universal (T-independent) 
part of the trapped virial coefficient is, 



b n ,oc (universal) = 



1 



,3/2 



Liiom 
n.oo " 



(12) 



Therefore, the higher density of states in traps suppresses 
the higher order virial coefficients, leading to an improved 
convergence of the expansion at low temperatures. 

High-T thermodynamics in free space. — Using rela- 
tion (fl2|) , we may determine the third virial coefficient in 
free space: b%°™ » -0.29095295. This does not agree 
with a previous field-theoretic calculation reported by 
Rupak [12], » 1.11. As well, we may calculate the 
high-T thermodynamics in free space, by taking a = 1 
and Qi = 2V{mk B T/2Tih 2 f/ 2 in Eqs. (pH and (ill- 
Fig. 4 presents the interaction energy Ei nt of a homo- 
geneous Fermi gas at unitarity as a function of temper- 
ature. For comparison, we also show the results of two 
quantum Monte-Carlo simulations. The virial expansion 
in free space seems to converge at T > 2Tp. 

In conclusion, we have shown that the virial expan- 
sion converges rapidly for a degenerate, resonant Fermi 
gas in a harmonic trap. This allows us to investigate 
the thermodynamics in a controllable way. We have pro- 
posed a practical method to obtain the third virial coeffi- 
cient throughout BCS-BEC crossover. Higher order coef- 
ficients are calculable in a similar manner, and may hold 



the prospect of revealing the exact thermodynamics 
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Figure 4: (color online). Temperature dependence of inter- 
action energy of a homogeneous gas at unitarity, obtained 
using different virial coefficients. For comparison, two quan- 
tum Monte Carlo data reported in Refs. @] and Q are shown. 
A finite range of interactions has been used in Ref. 0], which 
may lead to a systematic downshift in energies. We list also 
the result calculated by using Rupak's » 1.11. 



resonant Fermi gases in the deep degenerate regime. The 
current work provides a useful benchmark on testing fu- 
ture experiments and quantum Monte-Carlo simulations 
on strongly interacting Fermi gases. 
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